Pairing in exotic neutron rich nuclei around the drip line and in the crust of neutron 

stars 
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Exotic and drip-line nuclei as well as nuclei immersed in a low density gas of neutrons in the 
outer crust of neutron stars are systematically investigated with respect to their neutron pairing 
properties. This is done using Skyrme density-functional and different pairing forces such as a 
density-dependent contact interaction and a separable form of a finite-range Gogny interaction. 
Hartree-Fock-Bogoliubov and BCS theories are compared. It is found that neutron pairing is reduced 
towards the drip line while overcast by strong shell effects. Furthermore resonances in the continuum 
can have an important effect counterbalancing the tendency of reduction and leading to a persistence 
of pairing at the drip line. It is also shown that in these systems the difference between HFB and 
BCS approaches can be qualitatively large. 

PACS numbers: XXXXX 



I. INTRODUCTION 

Supcrfiuid Fermionic systems, in most of their realiza- 
tions in physics, are either infinite-size and uniform, like, 
e.g., superffuid 3, He or neutron matter in stars, or con- 
fined in a finite volume like, e.g., nuclei, cold atoms in 
traps, or metallic grains. In the former case, the pair- 
ing gap at the Fermi energy is a function of the density 
of matter while in the latter case, confinement induces a 
variation of the density on a scale which may be smaller 
or of the order of the coherence length of the Cooper 
pairs. For example the size of the Cooper pair in nuclei 
may vary locally by a big factor going from the size of 
the order of the nucleus in the interior to something like 
2 fm in the surface region [THl]- In cold atoms, the size 
of Cooper pairs can be varied with the help of Feshbach 
resonances and one can cover the whole range from the 
BCS to Bose-Einstein condensation (BEC) regimes. In 
nuclear physics, the local density approximation (LDA) 
is at its very limit of applicability because the size of 
the Cooper pairs is at the best, locally, of the size of the 
surface thickness Q. Nevertheless pairing correlations 
at the surface of finite nuclei may be expected to show 
some remembrance of what happens in infinite nuclear 
matter d, H, Q where the 1 So pairing gap is strongly 
peaked at a density close to po/5, where po is the sat- 
uration density of nuclear matter. Pairing correlations 
are, therefore, expected to be slightly enhanced at the 
surface of nuclei. In particular, one may naively think 
that in exotic or drip line nuclei neutron pairing is en- 
hanced in such situations. This would be due to the 
fact that the neutron density extends more smoothly out 
to low density and forms a more or less thick neutron 
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skin which may resemble a piece of low density neutron 
matter. This, however, is not the case as we will see in 
this work. The reason lies in the fact that the relation 
between infinite matter and surface of nuclei is, as just 
mentioned, not one-to-one, because of the long coherence 
length of the Cooper pairs As a consequence, finite- 
nuclei mostly reflect the pairing properties at a density 
that is the average density of the systems [Tf . The effect 
of the change in density on the pairing gap at the surface 
of nuclei is still difficult to analyze and to pin down from 
experiment [1,0. This, for instance, concerns the ques- 
tion whether and to what extent the pairing force may 
show a surface peaking flfj| . 

The main objective of the present work is, thus, the 
study of what happens to neutron pairing of nuclei 
around the neutron drip. This also concerns neutron 
rich nuclei at the frontier of stability, with a large neu- 
tron skin, as well as neutron rich nuclei embedded in a 
low density neutron gas as it can be in the inner crust 
of neutron stars. In this paper we are interested in what 
happens just at the overflow point where the neutron gas 
is about to appear, or has appeared at a very low density. 
Pairing correlations are built there in two rather different 
systems, the nuclear cluster and the shallow superfluid 
gas. Considered alone, these systems are paired, and put 
together, a mutual effect of pairing between the two com- 
ponents of the system might eventually change the pair- 
ing properties of the whole system. In cold atomic gases, 
it is possible to fabricate a trapping potential which goes 
over from a narrow container to a much wider one at 
a certain energy [III fl2l | . Filling up this potential with 
atoms one may also reach drip and overflow situations. 
Our study may, therefore, be of interest for other fields 
of physics as well. 

Since the limit of a nucleus embedded in a vanishing 
dilute gas is a drip-line nucleus, some properties of the 
latter type of nuclei, such as the existence of resonant 
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states in the continuum, can have an important impact 
on nuclei immersed in a low density gas (l3j . On the 
experimental side, nuclei at the border of stability cannot 
be created at the present time, but systematics towards 
neutron rich systems can be extracted from known nuclei 
masses. These systematics will therefore also be analyzed 
in the first part of the work. 

In this work an important issue is related to the ex- 
trapolation of the pairing properties of nuclear systems 
towards the limits of a very dilute external gas of neu- 
trons. This limit can be obtained either by increasing 
the size of the box at a fixed number of neutrons, or by 
varying the number of neutrons at a fixed size of the box. 
However, at the end, we will also study realistic Wigner 
Seitz cell scenarios in the context of the inner crust of 
neutron stars, where the neutron gas reaches non negli- 
gible densities. 

In addition, we also systematically investigate the men- 
tioned pairing properties using two different forces: a 
density-dependent contact interaction (DDCI) and a sep- 
arable finite-range interaction (SFRI) within the Hartree- 
Fock-Bogoliubov framework which is appropriate for in- 
homogeneous systems. The appropriateness of these two 
pairing interactions will be assessed. In inhomogeneous 
systems where the change in density is of the order of 
the coherence length or smaller, DDCI might be at its 
limit, since the interaction depends on the density via 
a local density approximation. It has been shown that, 
given a DDCI that reproduces the gaps of SFRI in both 
symmetric nuclear matter (SNM) and pure neutron mat- 
ter (PNM), the two interactions behave in a similar way 
also in inhomogeneous systems as nuclear clusters in the 
inner crust of neutron stars [l4| . 

In the present article, we discuss the signature of the 
superfluid state comparing two different theoretical pre- 
scriptions: we take the pairing gap, in canonical basis, 
closer to the Fermi energy, also called pairing gap of the 
Lowest Canonical State (LCS) Alcs 0, and the pair- 
ing gap Aw EH) averaged over all the states with the 
pairing tensor. We show that in finite systems like nu- 
clei or Wigner-Seitz cells, the pairing gap Alcs can be 
suppressed at overflow while the pairing energy and the 
pairing gap Ajjv may persists at overflow. The average 
pairing gap Ajjv and the pairing energy can, for certain 
superfluid features, be more appropriate quantities con- 
cerning the properties of inhomogeneous systems than 
the pairing gap A LC s- 

The paper is organized as follows: in Sec.[H]we present 
the equations we use to do our calculations and the 
methods to solve them; the results concerning nuclei 
around the drip-line are give in Sec. IIIIl while in Sec. IIVI 
we discuss the phenomenon of overflow in the passage 
outer/inner crust of a neutron star. Finally we give our 
conclusions in Sec. [Vl 



II. THE HARTREE-FOCK BOGOLIUBOV 
(HFB) THEORY 

The self-consistent HFB equations, see Eqs. @, are 



solved in a box on a spherical mesh with radius Rbox [12 1 ■ 
Using the standard notation (nlj, q) for the spherical 
single-particle states with radial quantum number n, or- 
bital angular momentum I, total angular momentum j, 
and isospin q—n, p, the single-particle wave functions 
(U, V) nlj ' q (r), are expanded on a basis of spherical Bessel 
functions, 



{U,V) nl ^(r) = J2(U,V)Z l ^u ail (r), 



(1) 



where u at i(r) — C a ^rji(k a jr), C a ^i is the normalization 
factor in the box, and ji are Bessel functions of the first 
kind with integer index I. The index a runs over a set of 
zeros of the Bessel function ji(k ai iRbox), going from the 
lower value k a =ij up to the momentum cutoff k max — 
4fm -1 . This corresponds to an HFB model-space energy 
cutoff of about h 2 k? nax /2m w 320 MeV (see Ref. El and 
references therein for more details). 

We use a Skyrme functional to build the single-particle 
Hamiltonian h and then we let the particles interact pair- 
wise in the pairing channel. The two-body matrix ele- 
ments of the pairing interaction in the J = 0, T = 1 
channel enter the neutron-neutron and proton-proton 
gap equations, whose solutions provide the matrix ele- 
ments of the state-dependent gap matrix A. The latter, 
in turn, enters the HFB equations, 
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(2) 

where fj, 9 F is the chemical potential and and V£ J ' q 

are the Bogoliubov amplitudes for the quasiparticle of 
energy E™ i< q . For the pairing channel we used two pair- 
ing interactions: a density dependent contact force and 
a finite range interaction in its separable approximation: 

(i) The two-body Density-Dependent Contact Interac- 
tion (DDCI) between particles at positions ri and 
r 2 reads El |IJ 



u(ri,r 2 ) = V 



Pb (R) 

Po 



S(r), 



(3) 



R = (ri + r 2 )/2 is the center of mass of the two 
interacting particles and r = ri — r 2 is their mu- 
tual distance. We choose V a = -530.0 MeV fm 3 , 
r] = 0.7, a = 0.45, p = 0.16 fm" 3 21]. We use 
a smooth cut-off acting in quasiparticle space at 
E a > 20 MeV, that is defined by an gaussian fac- 
tor exp(— (E a — 20) 2 /100), where a is a shorthand 
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notation for a — (nlj, q) . The parameters of this in- 
teraction are adjusted such that it mimics the gaps 
obtained using a Gogny force in SNM. 

(ii) A Separable Finite-Range pairing Interaction 
(SFRI) [H that reproduces the % Gogny D1S 
pairing gap at the Fermi surface in infinite nuclear 
matter (INM) [H, 



v(n, r 2 ,r' 1 ,4) = 7 P(r)P(r')<5(R - R')-(l - P° 



(4) 



The operator i(l — P a ) restricts the interac- 
tion to total spin S — 0. Strength and form 
factor are 7 = -728 MeVfrrT 3 and P(r) = 
l/(47r6 2 ) 3 / 2 exp(-r 2 /(46 2 )), where b = 0.644. The 
finite-range interaction in the pairing channel is 
added on top of a single-particle spectrum obtained 
with Skyrme interaction (see for instance Ref . [U] ) ■ 
A slight correction of the strength of the SFRI is 
therefore necessary and we fix 7 — > 0. 97 [23j . 

Since we are interested in systems at or beyond the 
neutron drip (i.e. a nucleus surrounded by a gas), we 
consider the same boundary conditions used in the cal- 
culation of Wigner Size cells in the inner crust of neutron 
stars. We thus impose the following Dirichlet-Neumann 
mixed boundary conditions Q : (i) even-parity wave func- 
tions vanish at R = R bo x] (ii) the first derivative of odd- 
parity wave functions vanishes at R = R box - When pre- 
senting the results for the HFB neutron pairing gap, we 
use two different definitions for the pairing gap. The 
first one, A^cSi is defined as the diagonal pairing matrix 
element corresponding to the canonical single-particle 
state (l5j . whose quasi- particle energy, 



(5) 



is the lowest. Here e a stands for the diagonal matrix 
element of the single-particle field /i' J ' 9 in canonical basis 
and A a the corresponding diagonal pairing-field matrix 
element. 

The second definition of the pairing gap, A^, is re- 
lated to the average of the state dependent gaps over the 
pairing tensor, i.e. 



A q 



(6) 



The pairing energy is defined in terms of the pairing 
tensor as, 

E q pair = \ £(2j + 1)£ U^A^V^. (7) 

nlj aa' 

A comparison of the pairing gaps A^cs an d Any with 
experimental value can be made considering some limi- 
tations due to other effects, such as the one induced by 



the time-odd term of the Hamiltonian. The experimen- 
tal gap for odd nuclei can, however, be deduced from the 
binding energies using a three-point formula centered on 
the odd nucleus, see the discussion in Ref. [25[ and ref- 
erences therein, as 

KZ( N ) = \ [E b (N + 1) - 2E b (N) + E b (N - 1)] , (8) 

where E b is the binding energy taken from Audi's 
database |2q |. The experimental pairing gap for even 
nuclei are deduced from the average of the three-point 
formula ((8|) applied to the two closest even nuclei as [13] , 

K7en(N) = \ - 1) + KZ(n + 1)] • (9) 

Definitions flSJ and (j9|) are closer to Alcs [28| than 
to Ajjv and, therefore, a comparison of the pairing gap 
A LC s an d the experimental gap © is shown in Fig. Q] 
for calcium, nickel, tin and lead isotopes. The error bars 
on the experimental gaps are estimated to be ±200 keV. 
This takes into account a small contribution coming from 
the experimental error bars on the masses, and a large 
contribution due to other contributions than the pairing 
gap on the experimental gap, such as for instance, the 
time odd terms in the mean field (29l - [3ll | , the use of the 3- 
point formula [3, H3| , 3-body terms in the pairing chan- 
nel [H, 33], and the particle- vibration coupling P, l33 . l35j . 

On the left side of Fig. [p a re shown the theoretical 
predictions based on SLy4 [24], [3^, 113] Skyrme interac- 
tion and SFRI in the pairing channel while on the ri ght , 
the results are obtained with the F+ interaction [38[ 
and SFRI. The evolution with the neutron number N of 
the experimental gap |9|) shown in Fig. [T] (top panels) is 
well reproduced by the HFB model with both the SLy4 
or F+ interaction in the mean-field and same interaction 
SFRI in the pairing channel. The theoretical calcula- 
tions have been performed up to the drip line, beyond 
the domain where experimental information are known. 
The limitations of the experimental data is quite visi- 
ble in Fig. [T] (bottom panels) where they are represented 
versus the asymmetry parameter (N — Z)/A. The ex- 
perimental data hardly reach (N — Z)/A w 0.25. It is 
interesting to notice that at the edge of the experimen- 
tal data the gaps (experimental and theoretical) tend to 
decrease with an important slope in the asymmetry di- 
rection (N — Z)/A, while beyond, the theoretical gaps 
go up again for larger (N — Z)/A > 0.3. The decrease 
of the experimental gaps ([9|) is mostly due to the shell 
closure at the boundary of experimental measurements. 
As a consequence, the asymmetry dependence of the ex- 
perimental pairing gap, which has been fitted as 



1 - 7.74 



N 



A 



6.75 

AV3 



(10) 



in the work of Ref. [39j, |40(, as well as in the former 
work of Ref. [4l|, tends to overestimate the asymmetry 
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FIG. 1: (Colors online) Graphical representation of experimental gaps for isotopes obtained with a 3-point filter and the 
theoretical gaps Alcs using the SLy4 (left panel) or F+ (right panel) Skyrme force and same separable pairing interaction. It 
should be underlined that the experimental error is obtained simply by the uncertainty of masses plus a constant value of 200 
KeV that is due to the limits of the theoretical model (see text for more details). In panels (c) and (d) the solid line corresponds 
to Eq. (|10[) . while the dashed line represents Eq. both lines been drawn for Sn isotopes. 



dependence of the pairing gap. This is clearly illustrated 
in Figs.QJc) andQJd) where the fit (fTCTj) is shown for Sn 
isotopes, i.e. Z=50 (solid line). Completing the unknown 
experimental data for large (N — Z) /A by the theoretical 
calculations shown in the upper panels of Fig. [TJ we can 
obtain a new set of parameters 



III. NUCLEAR SYSTEMS BEYOND THE 
DRIP-LINE 

In the following, the predictions for overflowing nuclear 
systems based on different pairing forces are analyzed. 



A = 



1 - 2 



N-Z 
A 



6.75 
~A^- 



(11) 



where the coefficient in front of the asymmetry is lower 
than in Eq. (ITU1) . The fit (ITT|) is also shown for Sn isotopes 
in Figs[IJc) and[IJd) (dashed lines). 



A. Global properties around the drip-line 

To investigate the behavior of different systems when 
crossing the drip line, we perform a systematic study 
of different isotopic chains, namely the Calcium, Nickel, 
Molybdenum, Tin, and Lead isotopes. The calculations 
have been done in a box of 40 fm radius and with the 
Dirichlet-Neumann mixed boundary conditions. 

In Fig. [5] we display the representative neutron gaps 
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FIG. 2: (Colors online) The evolution of pairing properties is shown along different isotopic chains: from top to bottom Calcium, 
Nickel and Molybdenum isotopes are represented. On the left are compared the two definitions of pairing gaps for Alcs and 
Auv using Skyrme SLy4 interaction for the ph channel and the SFRI ((4]) for the pairing channel. On the right, the pairing 
energy per neutron (0 is compared using different interactions in the ph channel such as F+, FO, F-, LNSI and SLy4. The 
vertical dashed line stands for the neutron drip-line nucleus using the same color as the associated interaction. 



Alcs and Ajjv (left panel) as well as the pairing en- 
ergy per neutron (right panel) corresponding to Calcium, 
Nickel, and Molybdenum isotopic chains computed us- 
ing the SFRI interaction for the pairin g ch annel and the 
mean-field provided by the SLy4, LNSI [IU, F+, F-, and 
FO [38[ Skyrme forces. In the left panel, the experimental 
neutron gaps are also displayed together with their asso- 
ciated error bars. The same quantities are shown in the 
two panels of Fig. |3jfor Tin and Lead isotopes. The iso- 
topes represented in Figs.f2]and[2]have been selected from 



their behavior at the drip line (43[: in Fig. [2] the isotopes 
are not magic at the drip line and some pairing correla- 
tions persist, while in Fig. [31 the isotopes are magic at the 
drip line and pairing correlations are strongly reduced at 
and beyond the drip. 

From Fig. [5] it can be observed that the two definitions 
of the pairing gap: Alcs [HI and Auv of Eq. ([6]) give 
quite similar predictions for bound nuclei within the ex- 
perimental error-bars, while they show a noticeable dif- 
ference at the drip-line and beyond. In particular it is 
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FIG. 3: (Colors online) As Fig[2] but for Tin and Lead isotopes. 



found that the gap Alcs is almost zero at the drip line 
for all nuclei analyzed in Figs.[2jand[3l while the gap Ajjv 
can persist with non zero values in some cases, namely for 
the Ca, Ni, and Mo isotopic chains shown in Fig. [51 while 
it is zero for other cases such as the Sn and Pb chains 
represented in Fig. [3j At the drip line, since the LCS-gap 
is the gap at the Fermi energy, the LCS-gap changes from 
its value in a bound state to the value corresponding to 
a delocalized unbound state at low density. It is known 
that the 1 So pairing gap in neutron matter is essentially 
zero at vanishing density Therefore, the LCS-gap is, 
indeed, expected to be quite suppressed at the drip-line. 

The pairing energies go to zero at the drip-line for Sn 
and Pb isotopes, see Fig. [3J while they persist in the case 
of Ca, Ni, and Mo as shown in Fig. [2] In Ref. [13], it has 
been argued that these differences are due to the presence 
of resonance states lying near the Fermi level of drip-line 
nuclei with non-negligible occupancy. In this case pairing 
correlations can persist and the pairing energy remains 
non-zero, as it happens in the Ca, Ni, and Mo isotopic 
chains. However, if there is a large energy gap between 
the last fully occupied bound state at the drip line and 
the first unoccupied resonant unbound state, as it is the 
case in the Sn and Pb isotopic chains, pairing correlations 
are strongly reduced at the drip-line. 

Since the pairing gaps Ajjy displayed on the left of 
Figs. [5J and [3J (see also Fig. [3} represent an average of 



the pairing correlations over all the states, it behaves 
similarly to the pairing energy see Eqs. ([6]) and (J7]). The 
qualitative difference between the pairing gaps Alcs and 
Ajjv shown in Fig. [5] is, therefore, simply related to the 
fact that in overflowing systems, there could be a signif- 
icant difference between the average pairing properties 
and the pairing gap corresponding to the last occupied 
state. The pairing gap which enters the ground state en- 
ergy is the pairing gap Ajjv, since it behaves like the 
pairing energy, while the pairing gap Alcs provides in- 
formation on the last occupied state which, e.g. influ- 
ences the gap in the level density and quantities which 
depend on it. The strong reduction of the pairing gap 
Alcs, being related to pairing property of a single state, 
does not necessary induce the suppression of the pairing 
energy, as shown in Fig. [51 Let us notice again that in 
stable nuclei, the pairing gaps Alcs and Ajjv are very 
similar (like in the nuclei represented on the left panels 
of Figs. [5J and [3]). It was, therefore, at first a surprise to 
observe a qualitative difference between the pairing gaps 
Alcs and Ajjy in overflowing systems. However, since 
Ajjv is an average of the gaps over the pairing tensor, 
it is clear that in regions where the individual Aj's vary 
rapidly, as it happens around the drip, an average will be 
different from the gap Alcs of the last occupied level. 
On the contrary, for stable nuclei, the individual gaps are 
smoothly varying and then Alcs ~ Auv- 
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Force 


m*/m 


m*/m 


Am* 


m*/m 


F+ 


0.700 


0.625 


0.170 


0.795 


FO 


0.700 


0.700 


0.001 


0.700 


F- 


0.700 


0.870 


-0.284 


0.586 


SLy4 


0.695 


0.800 


-0.186 


0.614 


LNS1 


0.604 


0.478 


0.342 


0.820 



TABLE I: In this table we show the isoscalar and isovector 
masses m*/m and m*/m, as well as the difference in neutron 
matter Am* = nCJm — rn^jm and the neutron effective mass 
ml/m for SLy4 0,113: F +, F ". p [H and LNS1 JH. 



In Figs. [2] and [3] are also shown the pai ring energies 
for several mean-field models, namely SLy4 [3 [S3, F+, 
F- 



F0 |38( and LNS1 |42j. It is observed that the reduc- 
tion of the gaps at the drip-line and beyond is a property 
which is independent of the considered models, while the 
absolute value of the pairing energy can vary from one 
model to another. The main difference among these mod- 
els is related to the effective mass in symmetric and in 
neutron matter. In this case the effective mass is defined 

H, 



as 



m 

TO* 



777, 

ml 



m 
m* 



rn 

rn* 



(12) 



where I — (p n — p v )j (p n + p p ) and q is the isospin charge 
{q = +1, —1 respectively for neutrons and protons). In 
Eq. (|12p the effective mass in symmetric matter is given 
by the isoscalar effective mass m*/m and the isovector 
effective mass to* /to is related to the isospin splitting 
in asymmetric matter. These quantities, as well as the 
effective mass splitting Am* = m^/m — m*/m and the 
neutron effective mass to* /to, both computed in neutron 
matter, are given in Tab. U for the Skyrme interactions 
represented in Figs. [2] and [3] 

For calcium isotopes, the pairing energy at the drip- 
line and beyond shown in Fig. [2] depends on the Skyrme 
interaction. There is indeed a group of Skyrme interac- 
tions for which the pairing gap is large at the drip line 
(F+, F0, F-), and another group for which it is reduced 
approximatively by a factor 2 at the drip line (SLy4, 
LNS1). This qualitative difference is due to the struc- 
ture of the drip line nuclei around the Fermi energy as 
it is shown in Tab.|n]for 66,68,70 Ca and for FQ ( top ) and 

SLy4 (bottom). The Skyrme interaction F0 predicts for 
66 Ca, the nucleus at the drip-line, a bound state <?g/2, 
in the canonical basis, close to the continuum with an 
energy -0.2 MeV, while SLy4 predict for 68 Ca, the 
nucleus at the drip line, a gg/2 state with lower energy 
—0.9 MeV. Since the g 9 / 2 is lower in energy with SLy4 
compared to F0, its occupation number is closer to 1 at 
the drip point and this state participates to a lower ex- 
tent to the pairing correlations. For both interaction, 



the system gains energy if instead of realizing the shell 
closure at N=50 and becoming non-superfluid, it leaves 
the g g /2 state partially unfilled so to gain extra pairing 
energy. This is what is observed in Fig. [2] not only at 
the drip line, but also beyond. The structure of the drip 
line nuclei is mostly conserved even in the presence of a 
gas. Despite these differences, it is interesting to remark 
that the global level occupation picture around the drip 
line is similar for the Skyrme interactions F0 and SLy4: 
the unbound states are occupied before the drip line is 
reached at variance with the usual claim that the drip 
occurs when the first unbound particle is produced [45| . 
In fact, among unbound states there are resonant states 
which play an important role around the drip line. They 
have large spatial overlaps with the single-particle lev- 
els of the nucleus and give a significant contribution to 
the pairing correlations. In most of the cases, resonant 
states are populated before unbound scattering states. 
They are very important to understand the transition 
from isolated nuclei to overflowing systems. 

We also analyzed the difference between the two kinds 
of pairing interactions which we used in this work: SFRI 
and DDCI. It is important to analyze the influence of 
the density dependence of the DDCI in inhomogeneous 
systems such as Wigner-Seitz cells at neutron overflow. 
A comparison of the pairing properties of isolated zirco- 
nium isotopes and in Wigner-Seitz cells is shown in Fig.0] 
for SFRI (left) and DDCI (right). The drip-line nucleus 
predicted by SLy4+SFRI is N = 84 while it is N = 88 for 
SLy4+DDCI. Despite this small difference, the behavior 
of the pairing gaps A^cs an d A^y, as well as the pair- 
ing energy is very similar at the drip line and beyond. 
We conclude that provided that the DDCI reproduce 
the same pairing gaps in symmetric and neutron mat- 
ter, DDCI give similar results compared to SFRI in in- 
homogeneous systems such as the Wigner-Seitz cells [i~4j . 
Otherwise the same scenario as in Figs. 2 and 3 is recov- 
ered also for the Zirconium isotopes. 



B. The limit of nuclei immersed in a vanishing 
dilute gas 

Below, we will explore what happens in the outer crust 
of neutron stars around the point where the neutrons 
start to drip into the free space between the lattice sites 
built by the nuclei. The question we want to answer is 
whether a very low density gas of supcrfluid neutrons in 
a large container can have any major influence on the 
superfluidity of the nuclei at the lattice sites. In this 
section, we, thus, explore, in a schematic study, the limit 
of nuclei immersed in a dilute gas with the density of the 
gas going to zero. 

We first show in Fig.[S]a zoom of Figs.[5]and[3]focussed 
on the overflowing nuclear systems. In calcium and nickel 
isotopes pairing correlations persists at overflow since the 
gap Auy do not vanish, while pairing correlations are 
almost suppressed at overflow in tin and lead isotopes. 
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TABLE II: Canonical single particle energies and canonical occupation probabilities for some calcium isotopes 




FIG. 4: (Colors online) Similar to Fig[2]but for Zr isotopes. On the first four panels the behavior using SLy4 Skyrme interaction 
and different pairing forces, DDCI (panel (a)-(b)) and SFRI (panel (c)-(d)), are compared. 
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FIG. 5: (Colors online) The gaps Alcs and Auv as a func- 
tion of the chemical potential. The calculations have been 
already shown in Fig. 12131 Here we make a zoom in the case 
of SLy4+SFRI for the four isotopic chains. On the same plot 
we show the calculation of Pure Neutron Matter (solid line). 



The solid line corresponds to the value of the pairing gap 
in uniform neutron matter for the densities of the gas. 
By construction, the gap A-lcs follows quite well the 
trend of the uniform neutron matter gap. However, in a 
realistic system, composed of a nucleus plus a gas, one 
has A[/y ^ A^cs an d, thus, Auv shows clear differences 
with the gas due to the influence of the resonance states 
in calcium and nickel isotopes. The understanding of 
overflowing systems, therefore, requires a better study of 
the pairing properties of nuclei immersed inside a gas at 
the limit of very low density. In the following we aim at 
decreasing the density of the gas to the lowest possible 
value. 

Starting from an overflowing system with an external 
gas, the low density limit can be reached in two different 
manners: the first one is by increasing the size of the 
box for a fixed number of neutrons and the second one 
is by decreasing the total number of neutrons at fixed 
box size. However, the numerical calculations cannot be 
performed in boxes with sizes larger than 80 fm with 
SFRI and 150 fm with DDCI. These limitations are due 
to the increasing number of partial waves as well as of the 
level density as the size of the box is increased. The larger 
size of the box reached with DDCI is related to the lower 
CPU time and memory request to perform calculations 
compared to SFRI. 

The effect of increasing the size of the box is illustrated 
in Fig. |B] for the case of 166 Zr. The total density and 
anomalous density profiles are represented with a linear 
and logarithmic scale. There is an important reduction of 
the gas density for boxes going from 20 to 100 fm, while 
the reduction of the density going from 100 to 150 fm is 
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FIG. 6: (Colors online) Evolution of the density p(r) (top 
panel) and abnormal density p(r) (bottom panel) for 166 Zr 
as a function of the box size for SLy4+SFRI. In the inset we 
show the semi-logarithmic scale. 



quite marginal. This is a limitation of this method which 
imposes to work with very large boxes to reach the low 
density limit. 

From the behavior of the density in 166 Zr as a function 
of the box radius represented in Fig. |H1 two regions can 
roughly be distinguished: one is the "bulk" and the other 
the gas. Fixing an arbitrary limit Ru m = 10 fm to sepa- 
rate the "bulk" from the "gas", the number of neutrons 
in the bulk can be estimated as, 



N 



bulk 



p„(r)d 3 r. 



(13) 



We obtain that for R box = 20 fm, 7V bulk w 99 neutrons 
and for R box = 100 fm, N hulk w 87. The number of neu- 
trons in the bulk decreases as a function of the box size 
having as a limit the isolated nucleus at the drip line. 
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FIG. 7: (Colors online) Evolution of the neutron pairing en- 
ergy for 166 Zr as a function of the box size for different func- 
tionals. The dashed lines represent the pairing energies of the 
last bound nucleus in the chain as shown in Tab. lIIII See text 
for details. 
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-1.47 
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-1.33 


0.000 
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TABLE III: In this table we show the exact value of the neu- 
tron chemical potential and pairing energy, expressed in MeV, 
for Zr isotopes shown in Fig|4] 



Since we perform constrained HFB calculation conserv- 
ing the total number of neutrons, the particles evapo- 
rated from the bulk appear in the scattering states. 

In Fig. [Jj we display by dots the difference be- 
tween the neutron pairing energy of the drip-line nu- 
cleus and that of the overflowing nuclear system 166 Zr, 
^pair( XZr ) ~ ^pair( 166Zr ): as a function of the box size 
for the models SLy4+DDCI, SLy4+SFRI and F0+DDCI. 
The considered drip-line nuclei in our models are 126 Zr for 
SLy4+DDCI, 124 Zr for SLy4+SFRI, 122 Zr for F0+DDCI. 
The values of the pairing energy in these nuclei are ex- 
tracted from Tab. IIIII The convergence to the asymp- 
totic value is different for the models SLy4+DDCI and 
SLy4+SFRI. We observe a fast convergence when going 
from a box of 20 fm to 40 fm for these two models. In such 
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FIG. 8: (Colors online) The pairing energy as a function of 
the external density of the neutron gas, Pg. The solid lines 
are obtained by fixing a box and increasing the number of 
particles from ±zo Zr to ^ uu Zr. The points are obtained by fix- 
ing the number of particles in 166 Zr and increasing the size of 
the box. The calculations are done using SLy4 and the DDDI 
pairing interaction 



case we have seen that the particle density is sufficiently 
high to completely fill the resonances and thus such states 
do not contribute to pairing superfluidity. The excess of 
pairing energy in those small boxes comes mainly from 
scattering states. In fact when we go from Rbox — 40fm 
to Rbox — 50 fm, the resonant state fi/ 2 starts to de- 
populate and thus we can form Cooper pairs using such 
state. This explains why we have an increase of super- 
fluidity. Going from Rb ox = 50 fm to Rb ox = 150 fm, the 
convergence becomes very slow. In this case the box is 
sufficiently large to decouple bound and scattering states, 
the residual pairing energy comes from the superfluidity 
of neutrons trapped into resonant states. The behavior 
is different in case of F0+DDCI model since a different 
mean field produces a different single particle structure. 

The asymptotic value is not reached in the larger boxes 
used in our calculations. Nevertheless, judging from the 
last numerical values, one may assume that the limit 
eventually goes to zero. 

We also explore the alternative scenario where the size 
of the box is kept fixed and the number of neutrons is 
varied from 88 to 160. The results computed for boxes 
of size ranging from 30 fm to 50 fm are displayed by 
solid lines in Fig. [5] These results should be compared 
with the ones shown in Fig. [7J To this end we show by 
points the pairing energy of the nucleus 166 Zr computed 
within boxes of different sizes as in Fig. [7j Again we ob- 
serve in Fig. 8 that for the case Rbox = 30 fm we have a 
quick drop of superfluidity at p™ ps 1.4 x 10 4 fm~ 3 and 
then an increase. This is the same phenomenon observed 
in Fig. [7j and it is due to a depopulation of the reso- 
nant state /V/2, that when occupied does not contribute 
to superfluidity and when is half filled gives an impor- 
tant contribution to superfluidity. In conclusion, the two 
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FIG. 9: (Colors online) We show the evolution of the canonical 
single particle levels, e n ij, of neutrons for the nucleus 166 Zr 
calculated using the SLy4 functional and a contact pairing 
force as a function of the box size. The dashed black line 
indicates the position of the Fermi energy. 

methods employed in this work to reach the low density 
limit in the gas show that the overflowing systems tends 
to the limit of the drip-line nucleus. The drip-line nu- 
cleus is therefore an important reference to understand 
and analyze the properties of the overflowing systems. 

In Fig. IH1 we show the evolution of the canonical sin- 
gle particle states for neutrons as a function of the box 
for 166 Zr. In this case, to simplify the picture, we de- 
cided to represent only the most important levels. When 
the size of the box increases, the level scheme is modi- 
fied. However, the resonance levels remain practically at 
the same position independently of the size of the box, 
as in the case of weakly bound nuclei [46| . With the 
help of this criterium we can identify the resonant state 
iV/2 located very close to the Fermi energy. Using the 
SLy4 plus DDCI interaction and performing the HFB 
calculation in a box of Rbox = 100 fm, we find that the 
occupancy of the / 7 / 2 resonant state is ~ 4.16 neutrons 
and the one of the P3/2 state is « 0.88 neutrons. From 
Tab. IIII1 we see that the corresponding last bound nu- 
cleus of the Zr-chain is 126 Zr. If we look at the canonical 
levels of this nucleus, we find a complete occupancy up 
to the state /111/2 at e^ 11/9 = —4.15 MeV, where we have 
a shell closure (N=82). There are two resonant states at 
e P3/2 = 0.653 MeV with occupation v^ 3/2 = 0.11 (corre- 
sponding to 0.44 neutrons) and e/ 7/2 = 0.15 MeV with 
occupation w^ 7/2 = 0.39 (corresponding to 2.73 neutrons). 
This means that we can find a stable nucleus although 
some of the particles stay in a resonant state close to 
zero. 

In Fig. 1101 we show the pairing energy per particle for 




FIG. 10: (Colors online) We compare the pairing energy per 
neutron calculated using the WS approximation for Mo iso- 
topes (circles) calculated as a function of the neutron Fermi 
energy ej. We observe that we can have both positive and 
negative values. For positive values we compare the results 
with a calculation of Pure Neutron Matter (PNM), solid line. 
The inset shows the behavior of the transition among bound 
and unbound nuclei. SLy4+SFRI. 



Mo isotopes as a function of the neutron chemical poten- 
tial /ip, using SLy4+SFRI model. The calculations have 
been performed for a fixed number of protons (Z=42) and 
fixed box radius (Rbox = 40 fm), similarly to what has 
been done by Grasso et al. [13] . We can observe in such a 
way the transition from bound nuclei to the gas+nucleus 
system. The results are compared with the analogous 
calculation done in neutron matter. We observe that for 
small positive values of fip the cluster+gas system has 
a bigger pairing energy than the homogeneous system. 
Such difference is due to the resonance structure as ex- 
plained in the previous sections. When we reach high 
density regions (fip > 2MeV), we see that the presence 
of the cluster reduces the pairing correlation. Such phe- 
nomenon has been already discussed by many authors 
concerning Wigner-Seitz calculations [3, [Tj, [U |4£h50| 



C. Detailed analysis of the resonant states 

In the previous sections, the role of resonant states 
have been stressed in order to understand the transition 
between nuclei and overflowing systems. To better de- 
scribe within a theoretical framework the resonant states 
we decided to solve the HFB equations in r-space treating 
in a proper way the continuum (without discretization). 

Defining, from the fully converged solution of Eq. , 
the hamiltonian h(R) and the pairing field A q (R) in the 
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FIG. 11: (Colors online) We show the occupation probability Vy(-E) for neutrons as a function of the quasi-particle energy for 
some given quantum numbers for 128 Zr (left) and 166 Zr (right) calculated using SLy4+DDCI functional starting from Eq. (|18p 
and for 4 given values of the size of the box in which we perform the calculations. See text for more details. 



following way, 




x ^(2j + l)U nl3 ' q {R)V nl3 ' q (R), (15) 



where W q {R) is the central potential and m*(R) is the 
effective mass. Eq. (TToT) is valid only in the case of a 
DDCI © while in the case of a finite range interaction, 
the pairing field is, in principle, a function of two vari- 
ables, see for instance Eqs. (IA4|) and (IB1[) . The pairing 
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field A q (R) entering into Eq. (|T5|) is, in such case, de- 
fined as A q (R) = A q jQC (R), where the local pairing field, 

A Loc( i ?)> is § iven b y 

Al OG (R)=A«(R,kUR)), (16) 

and the local Fermi momentum k F (R) is given by 

h 2 k 2 F {R) = 2m*(i?)(/i« - W*(R))). (17) 

where fi q F is the chemical potential also fixed by previous 
calculations. 

Eq. JD can now be expressed in coordinate space, fol- 
lowing Hamamoto et al. (5l| . 

(h(R) + t, F + E« p )U l i>«(R,E)-Ai(R)V l i>o(R,E) = 0, 
(h(R) + ^ F -E q p )V l ^ q (R,E)+A q (R)U^ q (R,E) = 0, 

(18) 

where the hamiltonian h(R) and the pairing field A q (R) 
are fixed from the converged solution of Eq. @ . 

In Eq. CD, the U lj > q (R, E) and V lj > q {R, E) amplitudes 
are the radial components of the quasi-particle wavcfunc- 
tions CI): see R- e f- and references therein for more 
details. Since only a single iteration is used to solve 
Eqs. CD, there is a small lack of consistency in this last 
calculation. It is, however, the price to pay for a proper 
treatment of the continuum states without discretization. 

We recall that in this case each quasi-particle energy E 
is an acceptable solution and, therefore, we can define an 
occupation probability as a function of the quasi-particle 
energy E as 

V 2 (E) = (2j + 1) J dRR 2 V l 2 j (E, R). (19) 

In Fig. [TT] we represent some of the solutions of Eq. CD 
for neutrons in 128 Zr (left panels) and 166 Zr (right pan- 
els) nuclei using the SLy4 Skyrme functional plus pairing 
with DDCI (JD • The label Rb ox in the upper left panel of 
Fig. [TT] stands for the size of the box in which the cen- 
tral potential W q (R) and the pairing field A q (R), which 
enter in Eq. CD, have been obtained. As expected, we 
observe that the energy of resonant states, such as p 3 / 2 , 
/ 5 / 2 , /V/2, in the drip line nucleus 128 Zr, panels (a), (c) 
and (e) of Fig. [TT] remain mostly constant as the size of 
the box Rbox increases. The situation can, however, be 
different in systems composed of a nucleus and a neu- 
tron gas, such as 166 Zr, since the density of the gas and 
the chemical potential are functions of Rbox- We notice, 
indeed, that the occupation of the resonant state P3/2 
for 166 Zr changes considerably when we change the size 
of the box, and this is due to the different position of 
the chemical potential compared to the position of the 
resonant state, see Fig. [§] The position of the centroid 
energy of the resonant state / 5 / 2 slightly increases as the 
size of the box increases. It reveals the sensitivity of the 
f 5 / 2 resonant state, being closed to the threshold energy, 
on the mean-field potential which is recalculated for each 



Rbox ■ It is however clear from the comparison of the left 
and right panels in Fig. [TTJ that the centroid energies of 
the resonant states p 3 / 2 , f^/2, fr/2 m 166 Zr (right panels) 
converge to the associated ones in 128 Zr (left panels) as 
the size of the box increases. 

We can as well define a total occupation probability in 
an energy region as 

V 2 = dEV^E). (20) 

Jo 

We display the value of this integral for different resonant 
states in 128 Zr and 166 Zr in the labels of Fig. [TT] These 
integrals have been calculated in an interval between 
and E max =5 MeV using different box sizes for W q (R) 
and A q (R). Comparing the results on the left and right 
sides of Fig. [TTJ we observe that there is a small fraction 
of particle that stays trapped in the resonances, that is 
the reason for the slow convergence of the evaporating gas 
to the value of the pairing gap of the last bound nucleus. 

From Fig. [TTJ we clearly see the resonant character of 
the state fr/2- Its centroid is located around «0.85 MeV 
and its width » 300 KeV for both 128 Zr and 166 Zr. The 
result is in good agreement with the canonical basis result 
shown in Fig. [9] Since this level is located very close to 
the Fermi energy we can use the approximate formula [5l| 

V(e - A) 2 + A 2 « A (21) 

The position of this peak mostly depends on the strength 
of the pairing force. It explains the stability of the reso- 
nant state ff/2 in 166 Zr as the size of the box increases. 

IV. PAIRING IN THE CRUST OF NEUTRON 
STARS 

The inner crust of neutron stars provides an excel- 
lent frame to ap ply the self-consistent mean-field the- 
ory 0, H3, Hj] . The inner crust extends from the drip 
density pdrip — 4 x 10 11 g cm~ 3 , where the neutrons start 
to leave from the nuclei into free space, till p~ 1.4x10 14 
g cm~ 3 , where the transition to uniform matter takes 
place. The inner crust of neutron stars is believed to be 
formed by a crystal lattice of nuclear clusters embedded 
in a low-density neutron gas and ultra-relativistic elec- 
trons. To describe crust matter, the Wigner-Seitz (WS) 
approximation is widely used. In this approximation the 
crust is divided into spherical cells, each one representing 
an inner crust region of a given average density. The WS 
cells are electrically neutral and the interaction among 
them is neglected in many cases. Since the seminal cal- 
culation of Negele and Vautherin in the inner crust of 
neutron stars [53| . more refined quantal calculations at 
HF or HFB level [1, [ij, of different degrees of 

complexity within the WS approximation have been per- 
formed. Also semiclassical models as the Constrained 
Liquid Drop Model [60( and Thomas- Fermi (TF) calcula- 
tions including pairing correlations [ill I6ll - |63T ] have been 
used to study the crust of neutrons stars. 
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TABLE IV: Equation of state of the outer crust based on Douchin-Haensel EoS [52| . 
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TABLE V: Equation of state of the inner crust based on Negele-Vautherin [53| ]. 



In the following, we will analyze the pairing properties 
of the lattice in the crust of neutron stars, based, for the 
outer crust, on the Douchin-Haensel equation of state 15211 
and, for the inner crust, on the Negele-Vautherin one [53j . 
The properties of the WS cells are listed in Table IIVI for 
the outer crust and Table |V] for the inner crust. We 
have not recalculated the WS configurations which min- 
imize the energy for each of the densities that we consid- 
ered and we have preferred, as a first step, to build the 
pairing correlations on WS configurations obtained from 
previous minimizations. Our choice is motivated by two 
reasons: first, we want to compare our results with other 
published previously, as in Refs. IHH, and second, 
we did not want to introduce self-consistency by varying 
at the same time the pairing interaction and the energy 
which would have introduced an important non-linear ef- 
fect in the search scheme. 

In Fig.[T21 we represent the local pairing gap A£ oc (i?), 



defined in Eq. (fTo) . for some WS cells representative 
of the inner crust, which are: 250 Zr, 500 Zr, 1100 Sn, and 
i800g n Q n .j- p p ane j f this figure we show the re- 
sults obtained by solving the full HFB equations given 
in Eq. ©, using the SLy4+SFRI model, while on the 
bottom panel we display the results obtained using the 
HF+BCS approximation, where only the diagonal cou- 
pling among pairing matrix elements in Eq.([2]) are con- 
sidered, by solid lines and the results computed with the 
TF+BCS approach [61] by dashed lines. For a presenta- 
tion of the Thomas-Fermi BCS approximation, see Ap- 
pendix [B] and references therein. We first discuss the 
HFB results shown on the top panel of Fig. [T2] The 
pairing field in the external region of the WS increases 
as the mass number of the WS cells increases. This is a 
well known phenomenon which can qualitatively be un- 
derstood from a local density approximation in the very 
low density regime: The density of the external neutron 
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FIG. 12: (Colors online) Local neutron pairing gap A£ oc (7?) 
for different WS cell calculated using the SLy4+SFRI. Lower 
panel: pairing gap as a function of neutron Fermi momentum, 
kp, for the PNM and the WS cells. See text for details. 



gas increases as the mass number of the cell goes up, 
and from uniform matter calculations, it is known that 
the average pairing gap increases as a function of the 
neutron density, for the cells considered in Fig. [121 We 
observe that in the gas region both, HFB, HF+BCS, and 
TF+BCS approaches give the same value for the pairing 
field Al oc (R). It is indeed expected that HFB and BCS 
theories coincide in uniform matter [l?], HU, [65| • 

The peak of the pairing field at the surface of the clus- 
ters can also be roughly justified from a LDA and it cor- 
responds to the maximum of the pairing gap in neutron 
matter (56| . However, the peak in LDA is quantitatively 
much higher than in the HF+BCS and TF+BCS calcu- 
lations |4ll66j|. 

The behavior of the pairing field inside the cluster is 
more complex to understand. From the HFB predictions, 
it is almost independent of the WS cells that we have 
considered, except for 250 Zr. In the case of 250 Zr, the 
reduction of the pairing field compared to the other cal- 
culations is induced by a shell effect: a small increase of 
the pairing strength makes the pairing field inside this 
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FIG. 13: (Colors online) Upper panel: pairing gap as a func- 
tion of neutron Fermi momentum, kp, for the PNM and the 
WS cells for the model SLy4+SFRI. Lower panel: the pairing 
energy per neutron. See text for details. 



cluster identical to that of the other WS cells. The in- 
dependence of the pairing field inside the cluster with 
respect to the outer gas density is typical of HFB the- 
ory as it can be seen in the upper panel of Fig. [T2] For 
the BCS approximations, shown on the bottom panel of 
Fig- Ell the pairing field inside the cluster increases as the 
mass number of the cell goes up. In the BCS approxima- 
tions, the pairing properties of the gas and of the cluster 
are strongly coupled, as observed in the bottom part of 
Fig. 1121 since pairing amplitudes U and V are diagonal in 
the HF basis. The non-diagonal pairing matrix elements 
in HFB theory strongly reduce the coupling between the 
gas and the cluster, as far as the pairing correlations are 
concerned. 

In Wigner-seitz cells, off-diagonal couplings in the pair- 
ing field play an important role in HFB theory, while they 
are neglected in the BCS approximation. Such a feature 
has been already remarked in finite nuclei [(37j where it 
was found that the use of the BCS approximation leads 
to a reduction of the gap compared to a full HFB solu- 
tion. We leave a better analysis concerning the difference 
between HFB and BCS to a future work. 
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In the upper panel of Fig. [13]we display the pairing gap 
at the Fermi energy for a realistic sequence of WS cells 
from the outer to the inner crust as shown in Tab. II VI and 
IVl using HFB (red circles) and TF+BCS approximation 
(blue squares), while in the lower panel we show the pair- 
ing energy per neutron in the same scenario. The pure 
neutron matter (PNM) solution of the BCS gap equation 
are also given in the two panels (solid line). HFB and 
TF+BCS approaches give similar pairing gaps and pair- 
ing energies in the region of positive chemical potential 
In the region of negative chemical potential fip, HFB 
gives larger pairing correlations than TF+BCS. This is 
basically due to the fact that the TF gap drops to zero 
at the neutron drip line [iH since in such semiclassical 
approaches shell effects are supressed. 

In the region of negative chemical potential fj,p, which 
coincides with the outer-crust, the contribution of the 
resonance states <i| and g| are given in Tab. I IVl It is 
shown that the resonance state d| plays an important 
role as the chemical potential approaches zero. A qual- 
itatively different behavior is observed around fip ~ 
between HFB and TF+BCS approaches, as illustrated 
in Fig. [T3J According to the TF+BCS approach, the 
pairing correlations tend to vanish around &s 0, as 
just mentioned and as previously observed in Ref. [111 ]. 
In the HFB theory, where resonance states as well as off- 
diagonal pairing matrix coupling are included, the pair- 
ing correlations around fXp » are smaller than at sta- 
bility but definitely not zero, as previously claimed in 
Ref. [43(. The fact that the TF approach gives a very 
much reduced gap at n n = may indicate a feature of 
the macroscopic limit. 

In the region of large positive chemical potentials, av- 
erage quantities such as A^cs an( i pairing energy are 
expected to be dominated by the neutron gas. However, 
differences in the local pairing gap inside the cluster (see 
Fig. n"2|) also induce some differences in the aforemen- 
tioned average quantities (see Fig. [T5)) . 

The effect emphasized in the previous sections, e.g. the 
coupling to resonant states close to the drip-line, is not 
clearly seen in Fig. Q2J The reason is that the change 
in neutron number is too sharp in the existing tables IIVI 
and El passing from 88 Ni (/i* < 0) to 180 Zr (^™ > 0). 
The pairing persistence phenomenon concerns only the 
first neutrons that start to drip out, as we have seen 
previously. It is now clear that a better investigation of 
the transition from the outer to the inner crust will be 
necessary in the future, using a smaller discretization on 
the values of the density of the equation of state. 



V. CONCLUSIONS 

In this work, we made a quite exhaustive study of pair- 
ing properties of nuclei around the neutron drip. In a 
LDA picture, one could have expected that neutron pair- 
ing is enhanced in such situations, since the neutron skin 
could be considered as a low density piece of neutron 



matter where pairing is at its maximum. However, nu- 
clei are too small and the pairing force too weak for LDA 
being a good description. The reality is more complex 
as our investigations show. Globally pairing is certainly 
not enhanced going from stability to the drip, rather it is 
reduced. However, the general feature is strongly hidden 
by nuclear shell effects. The isovector dependence of the 
pairing gap extrapolated from measured nuclei masses is 
too strong, and results as the consequence of an accumu- 
lation of closed shell nuclei at the border of the present 
experimental knowledge. From our theoretical calcula- 
tions, we predict that exotic neutron rich nuclei beyond 
this border shall exhibit a new raise of the pairing corre- 
lations compared to the present extrapolations yielding a 
much weaker average decrease to the drip than previously 
assumed [ioj . 

We should, moreover, mention that in this work our 
studies were restricted to spherical nuclei and that in 
this case and in most examples, the drip occurred at 
a magic, or close to magic neutron number with, nat- 
urally, reduced pairing correlations. The reason for 
this is not entirely clear. Either the nuclei search to 
gain binding energy in approaching magicity towards the 
drip because gain in energy by pairing is weakened, or 
it is the other way round, i.e., pairing is reduced be- 
cause anyway (spherical) nuclei drive to magicity at the 
drip. Furthermore, large scale nuclear mass calcula- 
tions, such as for instance Gogny D1S web mass table 
(http://www-phynu.cea.fr/HFB-Gogny.htm) show that 
only about half of the nuclei at the neutron drip are de- 
formed, a ratio much more in favor of sphericity than in 
the case of stable nuclei. In the deformed cases pairing 
acquires more usual values. Nevertheless, looking at the 
values given in Gogny D1S web mass table, pairing is 
certainly not enhanced with respect to the stable region 
in such cases either. These conclusions are based on the- 
oretical predictions but it can be surmised that reality is 
not entirely different. In any case, the situation of pairing 
properties of nuclei around the neutron drip is overcast 
by very large shell fluctuations as can be seen from the 
various figures given in the main text. An additional fea- 
ture which makes the situation complicated is the fact 
that there are resonances in the continuum which can be 
populated and which can have, in some cases, a sensible 
influence on the pairing properties of drip nuclei. This 
depends, for instance, on the precise position of the res- 
onances. In particular, for cases where a strongly degen- 
erate resonance level becomes located very closely to the 
chemical potential, pairing can become quite important 
like this is, e.g., the case for the Ca isotopes in Fig. [5J 

We have also analyzed the pairing correlations in the 
crust of neutron stars described in the Wigner-Seitz (WS) 
approach, which allowed us to study a scenario formed by 
nuclear clusters embedded in a low-density neutron gas. 
We see that in this situation the gap in the clusters is only 
weakly affected by the pairing in the neutron gas and that 
the cluster and the gas behave as almost independent 
systems. This result obtained from the state-of-the-art 
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HFB theory is not reproduced by its BCS approximation, 
since BCS is less adapted for such systems as WS cells as 
compared with stable nuclei, at least in the case of low 
neutron density. When the density of the gas decreases 
in approaching the transition to the outer crust where all 
the neutrons are bound, the average pairing gap and the 
pairing energy in the WS cell decrease following the trend 
of the neutron matter. The pairing correlations, though 
reduced at the transition between the inner-crust and 
the outer-crust, remain, however, non-zero. Our stud- 
ies have revealed that in this situation there is a large 
qualitative difference between HFB theory and semiclas- 
sical TF+BCS approximations. We postpone a detailed 
investigation of this feature to a future work. 
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Appendix A: Pairing field 



where u a .i(r) are defined in EclQ. 

As recently discussed in ref.[69| the 1 S component is 
by far the dominant one concerning calculations of pair- 
ing gaps at subnuclear densities. We can thus immedi- 
ately obtain the pairing field as 



A 9 (r*i, r 2 ) = -v{\n -r 2 \)$ q (fi,r 2 ) 



(A4) 



where v(\r\ — r 2 \) is the pairing interaction. Whose ma- 
trix elements can be written as 



(A5) 

To define a local pairing field, we need to apply the 
Wigner transformation [17j and write the pairing field 
as A q (R,k), where R is the center of mass coordinate 
and k is the relative momentum among two particles. 



Appendix B: Thomas-Fermi BCS approximation 

We present briefly the Thomas-Fermi BCS approxima- 
tion [6ll . |68[ . In such case we take the gap in phase space 
[riij at a momentum equal to the Fermi momentum as (65j 



A(R,k F ) = - 



dk' 



■v(k F -k')«;(R,k'). (Bl) 



Within the TF approach to the pairing problem, the 
above formula can still be written in a different way. 
Quantally the anomalous density matrix in r-space in the 
BCS approach is given by /t(r, r') = J2 n Kn ( r \ n )( n \ r ')- 
Therefore, after Wigner transformation, in the TF (h — > 
0) limit, one obtains 



In this appendix we briefly describe the numerical pro- 
cedure to obtain the pairing field in coordinate space. 
Following ref . [1] we can write it as 



$Vi,r 2 )= J2 



2 j + 1 



aa'nlj 



TTnlj,q\/nl3,q\ ,/,00 /•-> -f \ 
a V a> I Vaa'lj( r l^ r 2), 



(Al) 



where 4> a a >ij(ri,r 2 ) = [4>cxlj(ri)(f>a'ij(r 2 )]oo is the wave- 
function of two neutrons coupled to L = S = 0. 



^aa'ij(ri,r 2 ) = — 0„/j(ri)^ ;j (r 2 )F;( co sfi2)xoo, 

where xoo is the total spin function of two particles cou- 
pled to S=0. The wave function of the basis (fraijtfi) is 
defined as 



*»W= E <^^W^ S (A3) 



k(R,p) 



dEg 1 "{E)K{E)f E {^) 



(B2) 



where /e(R, p) is the normalized distribution func- 
tion [Z0|. Inserting Eq.(|B2} in Eq.(|Bl|), using the fact 
that in this case we can write the distribution function 
as /b(R, p) = S(E — H c i)/g(E) with H c i the classical 
Hamiltonian, and performing the angular average of the 
pairing force v(p,p') = J v(p — p')dQ (assuming that 
A(R, p) and /c(R, p) are spherically symmetric in mo- 
mentum space), the gap A(R, k = kp) can be finally re- 
cast as 



x J dEK(E)k E (R,)v(k F ,k E (R)), (B3) 

where k E (R) = ( 2m *i R) (E-V(R))^ 2 is the local Fermi 
momentum at energy E. 
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